Sección 4.7 - Ejercicio 10

a) Análisis descriptivo

El conjunto de datos Weekly contiene porcentajes de retorno semanales del índice S&P entre los años 1990 y 2010 con 1089 observaciones.

  • Year: es el año en el que fue obtenida la observación.
  • Lag 1-5: Retornos porcentuales de las 5 semanas anteriores, respectivamente.
  • Volume: Volumen de acciones intercambiadas (Número promedio de acciones diarias intercambiadas en billones).
  • Today: Retorno porcentual de la semana.
  • Direction: Indica si el mercado tuvo un retorno positivo o negativo en esta semana.

Encabezados y primeras 6 observaciones:

## # A tibble: 6 x 9
##    Year   Lag1   Lag2   Lag3   Lag4   Lag5 Volume  Today Direction
##   <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl>  <dbl> <fct>    
## 1  1990  0.816  1.57  -3.94  -0.229 -3.48   0.155 -0.27  Down     
## 2  1990 -0.27   0.816  1.57  -3.94  -0.229  0.149 -2.58  Down     
## 3  1990 -2.58  -0.27   0.816  1.57  -3.94   0.160  3.51  Up       
## 4  1990  3.51  -2.58  -0.27   0.816  1.57   0.162  0.712 Up       
## 5  1990  0.712  3.51  -2.58  -0.27   0.816  0.154  1.18  Up       
## 6  1990  1.18   0.712  3.51  -2.58  -0.27   0.154 -1.37  Down

Volumen promedio de acciones diarias intercambiadas en el tiempo

Se evidencia un crecimiento constante con un alza grande en el año 2005 y una estabilización de la trayectoria alrededor del año 2007. Se añadió un ruido a la posición de cada punto para poder evidenciar la densidad, además de transparencia. Se utilizó el método mgcv::gam para ajustar la tendencia.

El incremento en la volatilidad del índice acompañó su crecimiento a partir del 2007.

Correlación entre las variables

No se encontraron relaciones lineales claras entre las variables. Se presentan comportamientos esperados de la variable Today con la variable categórica Direction representada en el color, donde los valores negativos están asociados a un decrecimiento.

b) ¿Podemos predecir si el índice subirá o no mediante regresión logística?

Utilicemos regresión logística para predecir si habrá un retorno positivo del mercado en esta semana. Para ello, se utilizará el retorno porcentual de las 5 últimas semanas y el volumen sin separar un conjunto de entrenamiento y de prueba (Modelo sobreajustado).

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = weekly)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.6949  -1.2565   0.9913   1.0849   1.4579  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.26686    0.08593   3.106   0.0019 **
## Lag1        -0.04127    0.02641  -1.563   0.1181   
## Lag2         0.05844    0.02686   2.175   0.0296 * 
## Lag3        -0.01606    0.02666  -0.602   0.5469   
## Lag4        -0.02779    0.02646  -1.050   0.2937   
## Lag5        -0.01447    0.02638  -0.549   0.5833   
## Volume      -0.02274    0.03690  -0.616   0.5377   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1496.2  on 1088  degrees of freedom
## Residual deviance: 1486.4  on 1082  degrees of freedom
## AIC: 1500.4
## 
## Number of Fisher Scoring iterations: 4

Si realizamos una prueba de hipótesis formal para la significancia de los parámetros, encontraremos que la probabilidad de rechazar la hipótesis nula es muy grande para todas las variables predictoras excepto para el Lag2 donde al parecer hay significancia. El valor en el tiempo \(t-2\) parece tener una relación con el tiempo \(t\).

c) Matriz de confusión de la regresión logística y calidad de la predicción del modelo sobreajustado

##      Up
## Down  0
## Up    1

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary        0.561 
## 2 kap      binary        0.0350

Se obtiene una tasa de clasificación correcta accuracy del 56.1% utilizando regresión logística. Es importante interpretar este resultado correctamente, ya que la predicción obtenida indica que el mercado subió 557 días y bajó 54 días, y que no se está prediciendo el valor en el tiempo \(t+1\).

Cada predicción individual se debe interpretar de la siguiente manera: “Dadas las 5 semanas anteriores y el volumen de ésta, el mercado tendrá un retorno positivo”. La mayor cantidad de errores se cometió prediciendo que el mercado iba a subir y en realidad bajó, esta situación ocurrió 430 veces.

d) Modelo de regresión logística sólo con Lag2 y conjunto de prueba

## 
## Call:
## glm(formula = Direction ~ Lag2, family = "binomial", data = weekly_train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.536  -1.264   1.021   1.091   1.368  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  0.20326    0.06428   3.162  0.00157 **
## Lag2         0.05810    0.02870   2.024  0.04298 * 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1354.7  on 984  degrees of freedom
## Residual deviance: 1350.5  on 983  degrees of freedom
## AIC: 1354.5
## 
## Number of Fisher Scoring iterations: 4

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.625
## 2 kap      binary         0.141

Al igual que el modelo anterior, se tienen fallas al intentar predecir las bajas del mercado, sin embargo, este se validó en un conjunto de prueba de 104 observaciones desconocidas para el modelo. El accuracy en este caso es del 62.5%. Los resultados son sorpresivamente mejores que en el modelo sobreajustado del literal b).

d) Modelo LDA (Linear Discriminant Analysis) para predecir la tendencia del mercado

## Call:
## lda(Direction ~ Lag2, data = weekly_train)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581
## 
## Coefficients of linear discriminants:
##            LD1
## Lag2 0.4414162

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.625
## 2 kap      binary         0.141

Se obtienen las mismas clasificaciones para el modelo LDA y el modelo de regresión logística. Como el coeficiente discriminante lineal es Lag2 = 0.4414, es de esperar que en ambas categorías la gráfica de éstos sea similar. La matriz de confusión es igual a la del modelo logístico.

f) Modelo QDA (Quadratic Discriminant Analysis) para predecir la tendencia del mercado

## Call:
## qda(Direction ~ Lag2, data = weekly_train)
## 
## Prior probabilities of groups:
##      Down        Up 
## 0.4477157 0.5522843 
## 
## Group means:
##             Lag2
## Down -0.03568254
## Up    0.26036581

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.587
## 2 kap      binary         0

El clasificador QDA clasificó todos los valores como “Up”, si bien es una estrategia que otorga resultados buenos en términos porcentuales con respecto al accuracy, se espera que el modelo proponga en algunos casos predicciones de baja del mercado.

g) Modelo KNN (K-Nearest Neighbors) para predecir la tendencia del mercado

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary        0.548 
## 2 kap      binary        0.0453

El clasificador \(KNN\) tiene una tasa de clasificación correcta del 50% utilizando \(k = 1\), sin embargo, al utilizar \(k = 3\) se obtiene un 54.8%. Este clasificador no tuvo inconvenientes en proponer predicciones variadas según el valor anterior del mercado.

h) ¿Qué métodos producen los mejores resultados para predecir la tendencia del mercado utilizando Lag2?

Los métodos de regresión logística y análisis de discriminante lineal obtuvieron los mejores resultados, con un accuracy del 62.5%. Se decide favorecer al modelo de regresión logística dada la facilidad de interpretación de los coeficientes y las facilidades gráficas de comunicación que permite la curva logística.

i) Mejora del mejor modelo y variables adicionales de interés

Selección Step-wise con la función stepAIC sin interacciones

## 
## Call:  glm(formula = Direction ~ Lag1 + Lag2, family = "binomial", data = weekly_train)
## 
## Coefficients:
## (Intercept)         Lag1         Lag2  
##     0.21109     -0.05421      0.05384  
## 
## Degrees of Freedom: 984 Total (i.e. Null);  982 Residual
## Null Deviance:       1355 
## Residual Deviance: 1347  AIC: 1353

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary        0.577 
## 2 kap      binary        0.0350

El modelo seleccionado por Stepwise selection tiene un accuracy menor que el de un solo predictor, sin embargo, su AIC es el mismo. El AIC evalúa únicamente el ajuste, por lo que se recomendaría el modelo de una sola variable predictora. Es importante resaltar que el modelo de menor AIC podría ser mejor en un conjunto de prueba diferente.

Mejor modelo por tanteo de interacciones

Se proponen variables con interacciones y se selecciona el modelo con mayor accuracy en el conjunto de prueba, en este caso utilizando la variable Lag2 y la interacción entre Lag1, Lag 4 y Volume

## 
## Call:
## glm(formula = Direction ~ Lag2 + Lag1:Lag4:Volume, family = "binomial", 
##     data = weekly_train)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.493  -1.265   1.023   1.089   1.446  
## 
## Coefficients:
##                   Estimate Std. Error z value Pr(>|z|)   
## (Intercept)      0.2047966  0.0643540   3.182  0.00146 **
## Lag2             0.0553928  0.0291831   1.898  0.05768 . 
## Lag1:Lag4:Volume 0.0007316  0.0014747   0.496  0.61982   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1354.7  on 984  degrees of freedom
## Residual deviance: 1350.3  on 982  degrees of freedom
## AIC: 1356.3
## 
## Number of Fisher Scoring iterations: 4

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.644
## 2 kap      binary         0.179

Con base en los resultados anteriores se selecciona el modelo por tanteo por su accuracy superior, sin embargo, esto puede ser un sobreajuste del conjunto de prueba; por ello, si este modelo fuera a producción se recomienda realizar validación cruzada para identificar si la interacción de variables añadida es significativa. Si no es lo es, el modelo seleccionado por Stepwise selection se considera apropiado.

Sección 4.7 - Ejercicio 11

En este apartado, desarrollaremos un modelo para predecir cuando un carro tiene alto o bajo rendimiento de combustible respecto al kilometraje.

Inicialmente observemos las características de este conjunto de datos. El conjunto de datos Auto contiene información sobre 392 vehículos cuyas variables son:

  • mpg: millas por galón.
  • cylinders: numero de cilindros (entre 4 y 8)
  • displacement: desplazamiento del motor en pies.
  • horsepower: caballos de fuerza del motor.
  • weight: peso del vehículo en libras.
  • acceleration: tiempo que tarda en acelerar de 0 a 60 medido en segundos.
  • year: modelo del auto (módulo 100)
  • origin: origen del vehiculo.
  • name: nombre del vehículo.

Se mostrarán los primeros 6 datos junto con su encabezado.

## # A tibble: 6 x 9
##     mpg cylinders displacement horsepower weight acceleration  year origin
##   <dbl>     <dbl>        <dbl>      <dbl>  <dbl>        <dbl> <dbl>  <dbl>
## 1    18         8          307        130   3504         12      70      1
## 2    15         8          350        165   3693         11.5    70      1
## 3    18         8          318        150   3436         11      70      1
## 4    16         8          304        150   3433         12      70      1
## 5    17         8          302        140   3449         10.5    70      1
## 6    15         8          429        198   4341         10      70      1
## # ... with 1 more variable: name <fct>

a) Creación de una variable respuesta binaria

A continuación crearemos la variabel mpg01 que tomará el valor de 1 si la variable mpg esta por encima de la mediana, y 0 en caso contrario.

## # A tibble: 6 x 9
##   cylinders displacement horsepower weight acceleration  year origin name 
##       <dbl>        <dbl>      <dbl>  <dbl>        <dbl> <dbl> <fct>  <fct>
## 1         8          307        130   3504         12      70 1      chev~
## 2         8          350        165   3693         11.5    70 1      buic~
## 3         8          318        150   3436         11      70 1      plym~
## 4         8          304        150   3433         12      70 1      amc ~
## 5         8          302        140   3449         10.5    70 1      ford~
## 6         8          429        198   4341         10      70 1      ford~
## # ... with 1 more variable: mpg01 <fct>

b) Relaciones entre las variables

A continuación realicemos un análisis descriptivo para observar el comportamiento de las variables regresoras respecto a la variable mpg01

Por otro lado, demos una visualización para las variable origin, la cual es categórica.

Tomando en cuenta las gráficas, podemos afirmar que las variables cylinders, displacement, horsepower y weight son buenas variables cuantitativas para realizar una clasificación de la variable mpg01, pues al graficar discriminando por clases se puede observar que hay una particion notoria entre los carros de alto consumo y bajo consumo respecto a estas 4 variables. Note también que la variabele origen puede ser una buena variable predictora para mpg01 cuando el tipo de carro es americano.

c) División de conjuntos de entrenamiento y validación

A continuación dividiremos el conjunto Auto en dos conjuntos, uno de validación y otro de entrenamiento dejando el 75% de los datos para entrenar los modelos.

d) Modelo LDA

A continuación realizaremos un modelo LDA considerando las varables que fueron extraidas en el análisis descriptivo y que se consideraron importantes a la hora de realizar una clasificación para la variable respuesta mpg01.

## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight + 
##     origin, data = auto_train)
## 
## Prior probabilities of groups:
##         0         1 
## 0.5204082 0.4795918 
## 
## Group means:
##   cylinders displacement horsepower   weight    origin2    origin3
## 0  6.745098     270.0784  128.03922 3595.412 0.05882353 0.04575163
## 1  4.163121     115.5887   79.26241 2339.532 0.26241135 0.38297872
## 
## Coefficients of linear discriminants:
##                        LD1
## cylinders    -0.4317395853
## displacement -0.0018236085
## horsepower    0.0018307395
## weight       -0.0007360167
## origin2       0.3722004164
## origin3       0.5312690065

Grafiquemos a continuación el resultado del modelo.

Ahora construyamos la matriz de confusión para el modelo con el conjunto de datos de validación para observar el desempeño.

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.867
## 2 kap      binary         0.729

Note entonces que la tasa de clasificación correcta es de un 86.73%, lo que muestra que el modelo tiene una alta calidad.

e) Modelo QDA

Con el fin de realizar una comparativa con el modelo anterior, realicemos un modelo QDA para predecir la varibale mpg01 usando las variables más significativas encontradas en el análisis descriptivo.

## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower + weight + 
##     origin, data = auto_train)
## 
## Prior probabilities of groups:
##         0         1 
## 0.5204082 0.4795918 
## 
## Group means:
##   cylinders displacement horsepower   weight    origin2    origin3
## 0  6.745098     270.0784  128.03922 3595.412 0.05882353 0.04575163
## 1  4.163121     115.5887   79.26241 2339.532 0.26241135 0.38297872

Ahora construyamos la matriz de confusión para el modelo QDA con el conjunto de datos de validación.

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.867
## 2 kap      binary         0.729

Note entonces que la tasa de clasificación correcta es de un 86.73%, y en comparación con el modelo LDA, se puede decir que tienen un desempeño muy parecido de acuerdo a la metrica de la tasa de clasificación correcta.

f) Modelo de regresión logística

Por otra parte, construiremos un modelo de regresión logística para realizar una predicción de la variable mpg01.

## 
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower + 
##     weight + origin, family = binomial, data = auto_train)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -2.53196  -0.25902  -0.00542   0.33608   3.03362  
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  11.9583896  2.1576414   5.542 2.98e-08 ***
## cylinders    -0.1003094  0.4308485  -0.233  0.81590    
## displacement -0.0117784  0.0122590  -0.961  0.33665    
## horsepower   -0.0546728  0.0176520  -3.097  0.00195 ** 
## weight       -0.0015796  0.0009063  -1.743  0.08135 .  
## origin2       0.1646890  0.6822990   0.241  0.80927    
## origin3       0.6067756  0.6915181   0.877  0.38024    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 407.08  on 293  degrees of freedom
## Residual deviance: 152.08  on 287  degrees of freedom
## AIC: 166.08
## 
## Number of Fisher Scoring iterations: 7

Calculemos de igual manera la matriz de confusión para realizar una comparación con los métodos previos usando los datos de entrenamiento.

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.888
## 2 kap      binary         0.772

De acuerdo a este resultado, se puede ver que todos los métodos tienen un puntaje muy similar entre ellos pues en este caso, se puede observar una tasa de clasificación correcta del 88.78%.

g) Modelo \(KNN\)

Realicemos a continuación un modelo de \(KNN\) para predecir la variable mpg, en este caso, debemos probar con varios casos de \(K\) hallar un óptimo. Para esto, realizaremos una función que reciba como parámetro un numero \(k\) y esta función retornará la tasa de clasificación correcta para el \(k\) dado, realizando la clasificación sobre el conjunto de prueba.

Ahora tomaremos varios valores de \(k\) en un rango entre 1 y 60, y para estos valores, computaremos su tasa de clasificación correcta y los graficaremos para observar cual es el mejor valor de \(k\) en este rango.

Una vez observada la gráfica, encontraremos el mejor valor de \(k\) que mejora la tasa de clasificación correcta, y además para este valor de \(k\) realizaremos la matriz de confusión para ver el desempeño del modelo de manera global.

Concluimos entonces que para los valores de \(k\) probados, el mejor accuracy obtenido fue de 89.8%, el cual sigue siendo muy similar a los puntajes anteriores y desde un punto de vista global, estos puntajes denotan que los modelos tienen una muy buena calidad.

Sección 4.7 - Ejercicio 12

a) Implementación de la función Power()

A continuación implementaremos una función Power() que calculará la potencia \(2^3\)

b) Implementación de la función Power2() generalizada

Ahora, implementaremos la función Power2() que recibe dos parámetros x y a y computará \(x^a\)

c) Usos de la funcion Power2()

En el siguiente bloque de código calcularemos \(10^3\).

## [1] 1000

Ahora calcularemos el valor de \(8^{17}\).

## [1] 2.2518e+15

Y finalmente calcularemos el valor de \(131^3\).

## [1] 2248091

d) Implementación de la función Power3()

Ahora implementaremos una versión mucho más general para realizar las potencias que esta basada en la función Power2() pero en este caso, la función no imprimirá el resultado sino que retornara dicho valor como un objeto de R

e) Usando Plot3() para graficar la función \(f(x) = x ^ 2\)

En este bloque de código, graficaremos \(f(x) = x ^ 2\) usando una escala logaritmica para el eje \(y\)

f) Implementación de una forma generalizada para graficación

A continuación implementaremos un código para graficar la función \(f(x) = x ^ a\) tal que \(x \in I\) para un intervalo \(I \subseteq \mathbb{R}\).

Sección 4.7.1 - Ejercicio 13

Utilizaremos el conjunto de datos Boston para predecir si un suburbio tiene un porcentaje de crimen mayor o menor que la media. El conjunto de datos tiene las siguientes variables:

  • crim: Ratio de crimen per cápita por pueblo.
  • zn: Proporción de zonas residenciales asignadas para lotes mayores a 25.000 sq.ft.
  • indus: Proporción de negocios no comerciales en acres por pueblo.
  • chas: 1 si está en trayecto del río Charles y 0 en otro caso.
  • nox: Concentración de óxidos de nitrógeno (partes por millón).
  • rm: Número promedio de cuartos por hogar.
  • age: Proporción de unidades ocupadas construidas antes de 1940.
  • dis: Media ponderada de las distancias a 5 de los centros de empleo de Boston.
  • rad: Índice de accesibilidad a autopistas periféricas.
  • tax: Valor total de los impuestos por cada $10.000
  • ptratio: Ratio de pupilos-profesores por pueblo.
  • black: \(1000(Bk - 0.63)^2\) dónde \(Bk\) es la proporción de negros por pueblo.
  • lstat: Menor estatus de la población (porcentaje).
  • medv: Valor medio de una propiedad ocupada en $1000s.

Creación de la variable de interés

La mediana de la variable crim es 0.25651, se creará la variable crimen con niveles “alto” si es mayor a la mediana y “bajo” en otro caso.

Análisis descriptivo

Las variables de edad e indus presentan comportamientos discriminantes de interés, en la variable edad, los valores mayores a 70 tienen una gran prevalencia en la categoría alto y menores a este valor en la categoría bajo, por lo que al incluirla en uno de los modelos de clasificación se espera que sea significativa.

Este mismo comportamiento se identifica en la variable indus, donde los valores menores a 15 parecen tener una gran prevalencia en la categoría bajo, y entre 18 y 22 se encuentra la mayor cantidad de alto.

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Regresión logística

Construyamos inicialmente un clasificador únicamente con las dos características identificadas en el análisis descriptivo para evaluar las suposiciones. Realicemos también la separación en los conjutos de prueba y de entrenamiento.

## Joining, by = c("crim", "zn", "indus", "chas", "nox", "rm", "age", "dis", "rad", "tax", "ptratio", "black", "lstat", "medv", "crimen")
## 
## Call:
## glm(formula = crimen ~ age + indus, family = "binomial", data = boston_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.6774  -0.5407  -0.3414   0.5724   2.6981  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  4.727028   0.491097   9.625  < 2e-16 ***
## age         -0.047380   0.007062  -6.709 1.96e-11 ***
## indus       -0.131905   0.025355  -5.202 1.97e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 561.25  on 404  degrees of freedom
## Residual deviance: 337.58  on 402  degrees of freedom
## AIC: 343.58
## 
## Number of Fisher Scoring iterations: 5

Efectivamente ambos parámetros son altamente significativos. evaluemos la calidad de las predicciones.

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.772
## 2 kap      binary         0.547

Los resultados son aceptables, sin embargo, evaluemos el modelo step-wise propuesto para identificar otras variables de interés.

## 
## Call:
## glm(formula = crimen ~ . - crimen - crim, family = "binomial", 
##     data = boston_train)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -3.4328  -0.0032   0.0000   0.1377   1.6830  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  27.787319   8.586662   3.236  0.00121 ** 
## zn            0.062198   0.033850   1.837  0.06614 .  
## indus         0.068162   0.055431   1.230  0.21882    
## chas         -1.101626   0.829059  -1.329  0.18393    
## nox         -46.572894   8.258444  -5.639 1.71e-08 ***
## rm            0.079837   0.774442   0.103  0.91789    
## age          -0.033585   0.014415  -2.330  0.01981 *  
## dis          -0.659126   0.244784  -2.693  0.00709 ** 
## rad          -0.663217   0.167363  -3.963 7.41e-05 ***
## tax           0.006860   0.003033   2.262  0.02372 *  
## ptratio      -0.458391   0.151742  -3.021  0.00252 ** 
## black         0.035793   0.014010   2.555  0.01062 *  
## lstat        -0.006730   0.058173  -0.116  0.90790    
## medv         -0.141170   0.073683  -1.916  0.05538 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 561.25  on 404  degrees of freedom
## Residual deviance: 160.46  on 391  degrees of freedom
## AIC: 188.46
## 
## Number of Fisher Scoring iterations: 9
## 
## Call:  glm(formula = crimen ~ zn + nox + age + dis + rad + tax + ptratio + 
##     black + medv, family = "binomial", data = boston_train)
## 
## Coefficients:
## (Intercept)           zn          nox          age          dis  
##   26.279813     0.070380   -41.630025    -0.033514    -0.621363  
##         rad          tax      ptratio        black         medv  
##   -0.762505     0.008532    -0.402343     0.031772    -0.134618  
## 
## Degrees of Freedom: 404 Total (i.e. Null);  395 Residual
## Null Deviance:       561.2 
## Residual Deviance: 163.3     AIC: 183.3

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.901
## 2 kap      binary         0.800

El modelo stepwise tiene una mejora considerable en las métricas al incluir otras variables adicionales, sin embargo, se destaca que el poder discrimatorio de las variables adicionales con respecto al modelo inicial propuesto sólo mejoró la predicción en un 12.9%. En un ambiente de producción se recomendaría evaluar el costo de utilizar estas variables en el modelo versus el beneficio de predecir correctamente.

Modelo LDA

## Call:
## lda(crimen ~ zn + nox + age + dis + rad + tax + ptratio + black + 
##     medv, data = boston_train)
## 
## Prior probabilities of groups:
##      alto      bajo 
## 0.5111111 0.4888889 
## 
## Group means:
##             zn       nox      age      dis       rad      tax  ptratio
## alto  1.178744 0.6392367 86.25942 2.528803 14.787440 507.3575 18.97874
## bajo 22.209596 0.4703273 51.19040 5.149906  4.136364 306.2677 17.84192
##         black     medv
## alto 325.7682 20.33720
## bajo 390.5477 25.31162
## 
## Coefficients of linear discriminants:
##                  LD1
## zn       0.005331755
## nox     -8.025497066
## age     -0.016986785
## dis     -0.053038787
## rad     -0.074033448
## tax      0.001113865
## ptratio -0.062819427
## black    0.001227959
## medv    -0.033685626

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.891
## 2 kap      binary         0.778

Los resultados son muy similares a los del modelo de regresión logística siendo ligeramente peores, sin embargo, la gráfica construida por el LDA ofrece una interpretabilidad mayor a la obtenida en el literal 10.

Sección 8.4 - Ejercicio 7

En este literal, tomaremos el conjunto de datos Boston y analizaremos de una manera mas detallada el efecto de los parámetros en el \(MSE\).

Inicialmente creamos el conjunto de datos a analizar

## # A tibble: 6 x 14
##      crim    zn indus chas    nox    rm   age   dis   rad   tax ptratio
##     <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <int> <dbl>   <dbl>
## 1 0.00632    18  2.31 0     0.538  6.58  65.2  4.09     1   296    15.3
## 2 0.0273      0  7.07 0     0.469  6.42  78.9  4.97     2   242    17.8
## 3 0.0273      0  7.07 0     0.469  7.18  61.1  4.97     2   242    17.8
## 4 0.0324      0  2.18 0     0.458  7.00  45.8  6.06     3   222    18.7
## 5 0.0690      0  2.18 0     0.458  7.15  54.2  6.06     3   222    18.7
## 6 0.0298      0  2.18 0     0.458  6.43  58.7  6.06     3   222    18.7
## # ... with 3 more variables: black <dbl>, lstat <dbl>, medv <dbl>

Ahora graficaremos el error de prueba con los siguientes valores de mtry y ntree. Para calcular el error de prueba en cada uno de los valores, crearemos una función que evalua el error de prueba en cada caso.

Ahora dividiremos el conjunto de datos en datos de entrenamiento y datos de prueba. Se obtienen las siguientes dimensiones de los conjuntos de entrenamiento y prueba respectivamente.

## [1] 379  14
## [1] 127  14

Creamos un dataframe para los valores a probar

## # A tibble: 6 x 3
##   mtry_values ntree_values error
##         <int>        <int> <dbl>
## 1           2            3  26.5
## 2           3            3  23.5
## 3           4            3  16.9
## 4           5            3  21.2
## 5           6            3  16.9
## 6           7            3  15.2

Y finalmente graficamos los errores mediante una grafica de niveles como se muestra a continuación.

Observemos que para esta seleccion de conjuntos de entrenamiento y de validación, se cumple que el bagging no es el metodo mas efectivo, pues allí la grafica se muestra mas oscura, y por tanto evidencia un valor mas alto del MSE. Para este caso, el valor de mtree más adecuado se encuentra entre 4 y 5. Note que los valores de ntree no tienen una relacion alguna con el \(MSE\), pues los mejores valores del error de validación se pueden envontrar con valores de ntreecercanos a 100 pero tambien con valores cercanos a 50. Sin embargo, valores bajos del ntree dan errores muy altos sin importar el valor de mtree.

Sección 8.4 - Ejercicio 8

Analizaremos el conjunto de datos Carseats, el cuál es un conjunto de datos que posee información respecto a sillas para niño que se ubican dentro de los automoviles, y realizaremos una regresión para la variable Sales. Las variables que tiene este conjunto de datos son las siguientes:

  • Sales: unidades vendidas en miles en cada una de las ubicaciones.
  • CompPrice: precio por competidor en cada ubicación.
  • Income: ingresos de la comunidad en miles de dólares.
  • Advertising: presupuesto para propaganda local en cada ubicación.
  • Population: poblacion en cada región.
  • Price: precio de los asientos en cada sitio.
  • ShelveLoc: calidad de los puntos de venta en cada ubicación.
  • Age: edad promedio de la población en cada ubicación.
  • Education: nivel de educación en cada ubicación.
  • Urban: determina si la tienda se encuentra en territorio rural o urbano.
  • US: determina is la tienda esta en EE.UU.
## # A tibble: 6 x 11
##   Sales CompPrice Income Advertising Population Price ShelveLoc   Age
##   <dbl>     <dbl>  <dbl>       <dbl>      <dbl> <dbl> <fct>     <dbl>
## 1  9.5        138     73          11        276   120 Bad          42
## 2 11.2        111     48          16        260    83 Good         65
## 3 10.1        113     35          10        269    80 Medium       59
## 4  7.4        117    100           4        466    97 Medium       55
## 5  4.15       141     64           3        340   128 Bad          38
## 6 10.8        124    113          13        501    72 Bad          78
## # ... with 3 more variables: Education <dbl>, Urban <fct>, US <fct>

a) División del conjunto de datos en entrenamiento y validación

A continuación dividiremos el conjunto de datos en un conjunto de validación y otro de prueba, obteniendo las siguientes dimensiones respectivamente.

## [1] 300  11
## [1] 100  11

b) Entrenamiento de un árbol de regresión

Ahora, entrenaremos un árbol de regresión.

## node), split, n, deviance, yval
##       * denotes terminal node
## 
##  1) root 300 2383.00  7.412  
##    2) ShelveLoc: Bad,Medium 238 1373.00  6.706  
##      4) Price < 92.5 37  166.30  8.921  
##        8) Income < 47.5 7   37.60  6.526 *
##        9) Income > 47.5 30   79.18  9.480 *
##      5) Price > 92.5 201  992.00  6.298  
##       10) ShelveLoc: Bad 61  227.60  4.782  
##         20) Price < 129.5 43  145.30  5.262  
##           40) CompPrice < 125.5 26   64.85  4.477 *
##           41) CompPrice > 125.5 17   39.86  6.462 *
##         21) Price > 129.5 18   48.78  3.636 *
##       11) ShelveLoc: Medium 140  562.90  6.959  
##         22) Advertising < 14.5 118  406.30  6.626  
##           44) Price < 127 77  195.70  7.198  
##             88) CompPrice < 124.5 42   70.21  6.453 *
##             89) CompPrice > 124.5 35   74.12  8.093 *
##           45) Price > 127 41  138.20  5.553  
##             90) Age < 64.5 29   57.74  6.153 *
##             91) Age > 64.5 12   44.83  4.103 *
##         23) Advertising > 14.5 22   73.45  8.745  
##           46) Age < 54.5 14   25.93  9.592 *
##           47) Age > 54.5 8   19.86  7.261 *
##    3) ShelveLoc: Good 62  436.80 10.120  
##      6) Price < 109.5 19   59.60 12.420  
##       12) CompPrice < 120.5 13   17.09 11.570 *
##       13) CompPrice > 120.5 6   12.62 14.270 *
##      7) Price > 109.5 43  231.80  9.102  
##       14) Advertising < 13.5 35  150.70  8.500  
##         28) Age < 60.5 21   70.78  9.500 *
##         29) Age > 60.5 14   27.42  7.000 *
##       15) Advertising > 13.5 8   12.84 11.740 *

La gráfica del arbol se muestra a continuación.

Con este arbol entrenado, calculemos el error de prueba como se muestra a continuación

## [1] 4.477323

Y tambien mostremos el \(RMSE\) para este modelo tomando el conjunto de entrenamiento.

## [1] 2.115969

Como conclusión se puede observar que una variable que determina el numero de venta de manera diferenciadora es la variable ShelveLoc, note entonces que cuando un local tiene buena calidad en las estanterías del punto de venta, la cantidad de ventas es mucho mayor respecto a locales con una calidad mas baja. Aún así, los productos con un precio más bajo tiende a tener una cantidad de unidades vendidas mucho mas alta. Por tanto, dos aspectos que favorecen a las ventas es tener un local con buena calidad y precios por debajo de 97.5 dólares. Finalemente, considerando que la raiz del error cuadrático medio es aproximadamente de 2.1159687, quiere decir que el modelo se equivoca en un radio de 2.1159687, lo que realmente da una predicción aceptable sobre los asientos vendidos.

c) Uso de validacion cruzada para determinar el nivel de complejidad óptimo

En esta sección podaremos el árbol entrenado usando validación cruzada, con el objetivo de mejorar la calidad de la recomendación.

## $size
##  [1] 16 15 14 13 11 10  9  8  7  6  5  4  3  2  1
## 
## $dev
##  [1] 1498.228 1531.632 1529.217 1513.865 1516.933 1500.161 1484.006
##  [8] 1491.751 1547.555 1557.176 1556.349 1648.648 1873.875 1864.896
## [15] 2394.920
## 
## $k
##  [1]      -Inf  27.65922  29.89327  35.65000  37.04248  49.53619  51.34393
##  [8]  52.48000  68.26306  72.40865  83.19008 145.38605 201.45179 214.92784
## [15] 573.23081
## 
## $method
## [1] "deviance"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

Grafiquemos ahora la sumatoria de errores cuadráticos respecto al tamaño del árbol.

Note entonces que la validación cruzada arroja que el mejor árbol es el árbol original sin podar. Como se puede apreciar en la gráfica, el arbol mas grande es aquel que tiene mejor sumatoria de errores cuadráticos, por tanto no se recomienda podar el árbol.

d) Bagging e importancia de las variables

Ahora, usaremos bagging con el objetivo de observar si hay una mejora en la predicción y tambien se requiere observar cuáles son las variables mas importantes en este modelo.

Incialmente realizaremos bagging y compararemos el \(RMSE\) con los métodos previos.

## 
## Call:
##  randomForest(formula = Sales ~ ., data = carseats_train, importance = TRUE,      mtry = ncol(carseats) - 1) 
##                Type of random forest: regression
##                      Number of trees: 500
## No. of variables tried at each split: 10
## 
##           Mean of squared residuals: 2.551421
##                     % Var explained: 67.88

Calculemos ahora el error de validación para este método

## [1] 2.380399

Note entonces que el error de validación es mayor en este caso respecto al árbol de regresión. Ahora usemos el método importance() para determinar la importancia de las variables en este conjunto de datos.

##                %IncMSE IncNodePurity
## CompPrice   30.4360564    223.022071
## Income      11.1409987    147.056363
## Advertising 21.6230477    183.453069
## Population  -1.4934508     78.211847
## Price       62.9001945    676.637190
## ShelveLoc   71.8226799    754.037937
## Age         20.1073551    174.146019
## Education    5.5360020     71.905929
## Urban       -0.4547497      9.617858
## US           4.7873814     10.535014

Veamos una grafica de la importancia para cada una de las variables.

Note entonces que el incremento en la impureza y el porcentaje de incremento del MSE revelan que las variables ShelveLoc y Price son las que tienen mas relevancia dentro del modelo, y esta relevancia es bastante marcada respecto a las demás variables.

e) Uso de bósques aleatorios

En esta ocasión usaremos bosques aleatorios para determinar si este método es más efectivo y además, observaremos cuales son las variables de mayor importancia para este método considerando el conjunto de datos Carseats.

Primero entrenemos un modelo con los parámetros por defecto, obteniendo el siguiente \(MSE\)

## [1] 2.882621

Luego, analicemos la importancia de las variables usando la funcion importance()

##             IncNodePurity
## CompPrice       196.79119
## Income          198.82432
## Advertising     236.38564
## Population      147.20301
## Price           540.66745
## ShelveLoc       538.52677
## Age             234.85851
## Education       105.39935
## Urban            21.58035
## US               32.67248

Para tener una visualización de dicha importancia, consideremos la siguiente grafica.

Note entonces que, al igual que sucedio con el bagging, las dos variables mas importantes para este modelo son ShelveLoc y Price, y su importancia se resalta de manera significativa sobre las demás variables.

A continuación, veamos el efecto de \(m\) sobre el bosque aleatorio. Crearemos una función que reciba un valor de \(m\) y calcule el error de validación.

Ahora graficaremos el valor del \(MSE\) con respecto a \(m\).

Luego, extraeremos el mejor \(m\) como sigue.

## [1] 7

De esta manera, se muestra que considerando 7 variables dentro de la decisión en cada nodo, se obtiene un \(MSE\) de 2.394646.

Sección 8.4 - Ejercicio 9

El conjunto de datos OJ contiene 1070 compras donde los clientes compraron un jugo de naranja de Citrus Hill o de Maid Orange Juice y se guardaron un número de características de los clientes.

b) Árbol de decisión

## 
## Classification tree:
## tree(formula = Purchase ~ ., data = oj_train)
## Variables actually used in tree construction:
## [1] "LoyalCH"        "PriceDiff"      "WeekofPurchase" "ListPriceDiff" 
## Number of terminal nodes:  8 
## Residual mean deviance:  0.821 = 650.3 / 792 
## Misclassification error rate: 0.1925 = 154 / 800

El árbol resultante tiene 10 nodos terminales, la variable más importante es LoyalCH con una ponderación de 63, con una tasa de clasificación incorrecta de 14.88% en entrenamiento.

c) Interepretar el resultado de uno de los nodos

Interpretemos el nodo 1: En este nodo ingresan 800 observaciones del conjunto de entrenamiento, hay 311 observaciones de la clase CH y 489 de la clase MM, se realiza una partición a la izquierda con 512 observaciones y la derecha con 285 observaciones bajo el criterio \(LoyalCH < 0.48285\).

d) Gráfica del árbol y resultados

Los nodos utilizan las variables CH y PriceDiff para realizar la clasificación, los nodos terminales indican la clasificación final.

e) Predicciones, matriz de confusión y clasificación correcta

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.830
## 2 kap      binary         0.649

Se obtiene una tasa de clasificación correcta del 80%. En general la mayor cantidad de errores se cometió al clasificar las variables CH como MM, esto ocurrió 29 veces.

f) Aplicar cv.tree() para determinar el tamaño óptimo del árbol

## $size
## [1] 8 7 4 2 1
## 
## $dev
## [1] 188 187 187 189 308
## 
## $k
## [1]        -Inf   0.0000000   0.6666667   2.5000000 146.0000000
## 
## $method
## [1] "misclass"
## 
## attr(,"class")
## [1] "prune"         "tree.sequence"

g) Gráfica con el tamaño del árbol y la validación cruzada

h) Árbol con la menor tasa de error de validación cruzada incorreta

De acuerdo a la gráfica anterior, el árbol con tamaño 5 parece obtener el menor error de validación cruzada. Realicemos la poda de acuerdo a esta inforamción y grafiquemos el nuevo árbol.

j) Comparación de los errores de clasificación en entrenamiento

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.808
## 2 kap      binary         0.598

Se obtiene una clasificación correcta del 84.5% para el modelo CART podado, esto se compara con el modelo completo que da una clasificación correcta del 85.2%, es decir, el modelo con validación cruzada es ligeramente peor en el conjunto de entrenamiento pero tiene una mayor capacidad de generalización gracias a que evita sobreajustar los datos.

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.808
## 2 kap      binary         0.598

k) Comparación los resultados del árbol podado vs el completo en validación

## # A tibble: 2 x 3
##   .metric  .estimator .estimate
##   <chr>    <chr>          <dbl>
## 1 accuracy binary         0.830
## 2 kap      binary         0.649

En el árbol sin podar obtuvimos una tasa de clasificación correcta del 80%, sin embargo, al realizar la poda se obtuvo una tasa de clasificación correcta del 80.4%, que podría ser una mejora apreciable dependiendo del volumen de ventas de jugo de naranja asociado que trae la predicción correcta.

Sección 8.4 - Ejercicio 10

El objetivo de esta sección es usar boosting para hacer una predicción sobre la variable Salary en el conjunto de datos Hitters

El conjunto de datos recoge información sobre jugadores de béisbol de las ligas mayores de Estados Unidos. A continuación se muestra un resumen de dicho conjunto de datos.

## # A tibble: 6 x 20
##   AtBat  Hits HmRun  Runs   RBI Walks Years CAtBat CHits CHmRun CRuns  CRBI
##   <int> <int> <int> <int> <int> <int> <int>  <int> <int>  <int> <int> <int>
## 1   293    66     1    30    29    14     1    293    66      1    30    29
## 2   315    81     7    24    38    39    14   3449   835     69   321   414
## 3   479   130    18    66    72    76     3   1624   457     63   224   266
## 4   496   141    20    65    78    37    11   5628  1575    225   828   838
## 5   321    87    10    39    42    30     2    396   101     12    48    46
## 6   594   169     4    74    51    35    11   4408  1133     19   501   336
## # ... with 8 more variables: CWalks <int>, League <fct>, Division <fct>,
## #   PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>,
## #   NewLeague <fct>

a) Remoción de NA’s y escalamiento logaritmico

Inicialmente se eliminaran aquellas filas donde el salario tenga el valor NA. Para ello primero observemos si la columna Salary es la única que tiene NA.

##      AtBat            Hits         HmRun            Runs       
##  Min.   : 16.0   Min.   :  1   Min.   : 0.00   Min.   :  0.00  
##  1st Qu.:255.2   1st Qu.: 64   1st Qu.: 4.00   1st Qu.: 30.25  
##  Median :379.5   Median : 96   Median : 8.00   Median : 48.00  
##  Mean   :380.9   Mean   :101   Mean   :10.77   Mean   : 50.91  
##  3rd Qu.:512.0   3rd Qu.:137   3rd Qu.:16.00   3rd Qu.: 69.00  
##  Max.   :687.0   Max.   :238   Max.   :40.00   Max.   :130.00  
##                                                                
##       RBI             Walks            Years            CAtBat       
##  Min.   :  0.00   Min.   :  0.00   Min.   : 1.000   Min.   :   19.0  
##  1st Qu.: 28.00   1st Qu.: 22.00   1st Qu.: 4.000   1st Qu.:  816.8  
##  Median : 44.00   Median : 35.00   Median : 6.000   Median : 1928.0  
##  Mean   : 48.03   Mean   : 38.74   Mean   : 7.444   Mean   : 2648.7  
##  3rd Qu.: 64.75   3rd Qu.: 53.00   3rd Qu.:11.000   3rd Qu.: 3924.2  
##  Max.   :121.00   Max.   :105.00   Max.   :24.000   Max.   :14053.0  
##                                                                      
##      CHits            CHmRun           CRuns             CRBI        
##  Min.   :   4.0   Min.   :  0.00   Min.   :   1.0   Min.   :   0.00  
##  1st Qu.: 209.0   1st Qu.: 14.00   1st Qu.: 100.2   1st Qu.:  88.75  
##  Median : 508.0   Median : 37.50   Median : 247.0   Median : 220.50  
##  Mean   : 717.6   Mean   : 69.49   Mean   : 358.8   Mean   : 330.12  
##  3rd Qu.:1059.2   3rd Qu.: 90.00   3rd Qu.: 526.2   3rd Qu.: 426.25  
##  Max.   :4256.0   Max.   :548.00   Max.   :2165.0   Max.   :1659.00  
##                                                                      
##      CWalks        League  Division    PutOuts          Assists     
##  Min.   :   0.00   A:175   E:157    Min.   :   0.0   Min.   :  0.0  
##  1st Qu.:  67.25   N:147   W:165    1st Qu.: 109.2   1st Qu.:  7.0  
##  Median : 170.50                    Median : 212.0   Median : 39.5  
##  Mean   : 260.24                    Mean   : 288.9   Mean   :106.9  
##  3rd Qu.: 339.25                    3rd Qu.: 325.0   3rd Qu.:166.0  
##  Max.   :1566.00                    Max.   :1378.0   Max.   :492.0  
##                                                                     
##      Errors          Salary       NewLeague
##  Min.   : 0.00   Min.   :  67.5   A:176    
##  1st Qu.: 3.00   1st Qu.: 190.0   N:146    
##  Median : 6.00   Median : 425.0            
##  Mean   : 8.04   Mean   : 535.9            
##  3rd Qu.:11.00   3rd Qu.: 750.0            
##  Max.   :32.00   Max.   :2460.0            
##                  NA's   :59

Note entonces que la única columna que tiene NA es Salary. Ahora podemos eliminar estas filas y verificar que no quedan NA

##      AtBat            Hits           HmRun            Runs       
##  Min.   : 19.0   Min.   :  1.0   Min.   : 0.00   Min.   :  0.00  
##  1st Qu.:282.5   1st Qu.: 71.5   1st Qu.: 5.00   1st Qu.: 33.50  
##  Median :413.0   Median :103.0   Median : 9.00   Median : 52.00  
##  Mean   :403.6   Mean   :107.8   Mean   :11.62   Mean   : 54.75  
##  3rd Qu.:526.0   3rd Qu.:141.5   3rd Qu.:18.00   3rd Qu.: 73.00  
##  Max.   :687.0   Max.   :238.0   Max.   :40.00   Max.   :130.00  
##       RBI             Walks            Years            CAtBat       
##  Min.   :  0.00   Min.   :  0.00   Min.   : 1.000   Min.   :   19.0  
##  1st Qu.: 30.00   1st Qu.: 23.00   1st Qu.: 4.000   1st Qu.:  842.5  
##  Median : 47.00   Median : 37.00   Median : 6.000   Median : 1931.0  
##  Mean   : 51.49   Mean   : 41.11   Mean   : 7.312   Mean   : 2657.5  
##  3rd Qu.: 71.00   3rd Qu.: 57.00   3rd Qu.:10.000   3rd Qu.: 3890.5  
##  Max.   :121.00   Max.   :105.00   Max.   :24.000   Max.   :14053.0  
##      CHits            CHmRun           CRuns             CRBI       
##  Min.   :   4.0   Min.   :  0.00   Min.   :   2.0   Min.   :   3.0  
##  1st Qu.: 212.0   1st Qu.: 15.00   1st Qu.: 105.5   1st Qu.:  95.0  
##  Median : 516.0   Median : 40.00   Median : 250.0   Median : 230.0  
##  Mean   : 722.2   Mean   : 69.24   Mean   : 361.2   Mean   : 330.4  
##  3rd Qu.:1054.0   3rd Qu.: 92.50   3rd Qu.: 497.5   3rd Qu.: 424.5  
##  Max.   :4256.0   Max.   :548.00   Max.   :2165.0   Max.   :1659.0  
##      CWalks       League  Division    PutOuts          Assists     
##  Min.   :   1.0   A:139   E:129    Min.   :   0.0   Min.   :  0.0  
##  1st Qu.:  71.0   N:124   W:134    1st Qu.: 113.5   1st Qu.:  8.0  
##  Median : 174.0                    Median : 224.0   Median : 45.0  
##  Mean   : 260.3                    Mean   : 290.7   Mean   :118.8  
##  3rd Qu.: 328.5                    3rd Qu.: 322.5   3rd Qu.:192.0  
##  Max.   :1566.0                    Max.   :1377.0   Max.   :492.0  
##      Errors           Salary       NewLeague
##  Min.   : 0.000   Min.   :  67.5   A:141    
##  1st Qu.: 3.000   1st Qu.: 190.0   N:122    
##  Median : 7.000   Median : 425.0            
##  Mean   : 8.593   Mean   : 535.9            
##  3rd Qu.:13.000   3rd Qu.: 750.0            
##  Max.   :32.000   Max.   :2460.0

Ahora realicemos un escalamiento logaritmico sobre la columna Salary

## # A tibble: 6 x 20
##   AtBat  Hits HmRun  Runs   RBI Walks Years CAtBat CHits CHmRun CRuns  CRBI
##   <int> <int> <int> <int> <int> <int> <int>  <int> <int>  <int> <int> <int>
## 1   315    81     7    24    38    39    14   3449   835     69   321   414
## 2   479   130    18    66    72    76     3   1624   457     63   224   266
## 3   496   141    20    65    78    37    11   5628  1575    225   828   838
## 4   321    87    10    39    42    30     2    396   101     12    48    46
## 5   594   169     4    74    51    35    11   4408  1133     19   501   336
## 6   185    37     1    23     8    21     2    214    42      1    30     9
## # ... with 8 more variables: CWalks <int>, League <fct>, Division <fct>,
## #   PutOuts <int>, Assists <int>, Errors <int>, Salary <dbl>,
## #   NewLeague <fct>

b) Seleccion de los conjuntos de prueba y validación

Se seleccionan las primeras 200 filas como el conjunto de prueba y las filas restantes se seleccionan como conjunto de validación, obteniendo los dos siguientes conjuntos de datos de manera respectiva.

## # A tibble: 200 x 20
##    AtBat  Hits HmRun  Runs   RBI Walks Years CAtBat CHits CHmRun CRuns
##    <int> <int> <int> <int> <int> <int> <int>  <int> <int>  <int> <int>
##  1   315    81     7    24    38    39    14   3449   835     69   321
##  2   479   130    18    66    72    76     3   1624   457     63   224
##  3   496   141    20    65    78    37    11   5628  1575    225   828
##  4   321    87    10    39    42    30     2    396   101     12    48
##  5   594   169     4    74    51    35    11   4408  1133     19   501
##  6   185    37     1    23     8    21     2    214    42      1    30
##  7   298    73     0    24    24     7     3    509   108      0    41
##  8   323    81     6    26    32     8     2    341    86      6    32
##  9   401    92    17    49    66    65    13   5206  1332    253   784
## 10   574   159    21   107    75    59    10   4631  1300     90   702
## # ... with 190 more rows, and 9 more variables: CRBI <int>, CWalks <int>,
## #   League <fct>, Division <fct>, PutOuts <int>, Assists <int>,
## #   Errors <int>, Salary <dbl>, NewLeague <fct>
## # A tibble: 63 x 20
##    AtBat  Hits HmRun  Runs   RBI Walks Years CAtBat CHits CHmRun CRuns
##    <int> <int> <int> <int> <int> <int> <int>  <int> <int>  <int> <int>
##  1   419   101    18    65    58    92    20   9528  2510    548  1509
##  2   376    82    21    42    60    35     5   1770   408    115   238
##  3   486   145    11    51    76    40    11   3967  1102     67   410
##  4   246    76     5    35    39    13     6    912   234     12   102
##  5   205    52     8    31    27    17    12   5134  1323     56   643
##  6   348    90    11    50    45    43    10   2288   614     43   295
##  7   523   135     8    52    44    52     9   3368   895     39   377
##  8   312    68     2    32    22    24     1    312    68      2    32
##  9   496   119     8    57    33    21     7   3358   882     36   365
## 10   126    27     3     8    10     5     4    239    49      3    16
## # ... with 53 more rows, and 9 more variables: CRBI <int>, CWalks <int>,
## #   League <fct>, Division <fct>, PutOuts <int>, Assists <int>,
## #   Errors <int>, Salary <dbl>, NewLeague <fct>

c) Aplicación de boosting y determinación del parametro \(\lambda\) de contracción

En esta sección realizaremos boosting y calcularemos el mejor parámetro de contracción \(\lambda\) sobre un determinado rango. Para esto, crearemos una función que calcule el \(MSE\) recibiendo \(\lambda\) como parámetro.

Ahora definamos el rango sobre el cual se quiere graficar los \(MSE\), y realicemos dicha gráfica usando la función que se definió anteriormente.

En esta gráfica se puede apreciar que el mejor valor para \(\lambda\) es 0.17 con el cual se obtiene un \(MSE\) de 0.0446519.

d) Comparación de boosting con métodos previos

En esta sección compararemos el \(MSE\) de boosting respecto a otros métodos vistos ateriormente como regresión lineal múltiple y regresión de componentes principales.

Primero calculemos el \(MSE\) para el metodo de regresión lineal múltiple

## [1] 0.0188556

Por otro lado, apliquemos regresión con componentes principales para extraer el \(MSE\).

## Data:    X dimension: 200 19 
##  Y dimension: 200 1
## Fit method: svdpc
## Number of components considered: 19
## 
## VALIDATION: RMSEP
## Cross-validated using 10 random segments.
##        (Intercept)  1 comps  2 comps  3 comps  4 comps  5 comps  6 comps
## CV          0.3981   0.2781   0.2806   0.2798   0.2799   0.2786   0.2761
## adjCV       0.3981   0.2778   0.2802   0.2792   0.2794   0.2779   0.2755
##        7 comps  8 comps  9 comps  10 comps  11 comps  12 comps  13 comps
## CV      0.2756   0.2738   0.2741    0.2758    0.2777    0.2774    0.2780
## adjCV   0.2750   0.2732   0.2735    0.2750    0.2769    0.2765    0.2768
##        14 comps  15 comps  16 comps  17 comps  18 comps  19 comps
## CV       0.2742    0.2768    0.2710    0.2706    0.2733    0.2749
## adjCV    0.2731    0.2756    0.2699    0.2693    0.2718    0.2733
## 
## TRAINING: % variance explained
##         1 comps  2 comps  3 comps  4 comps  5 comps  6 comps  7 comps
## X         39.40    60.67    70.76    79.64    84.73    89.29    92.38
## Salary    52.42    52.67    53.22    53.56    54.56    55.34    55.70
##         8 comps  9 comps  10 comps  11 comps  12 comps  13 comps  14 comps
## X         95.02    96.32     97.26     98.06     98.72     99.23     99.54
## Salary    56.48    56.68     56.77     56.80     57.41     57.91     58.68
##         15 comps  16 comps  17 comps  18 comps  19 comps
## X          99.77     99.90     99.97     99.99    100.00
## Salary     58.75     60.27     61.27     61.27     61.49

En este caso, encontramos que el menor valor optimo se logró con \(M = 17\) por tanto calcularemos el error de validación con dicho número de componentes, obteniendo así el siguiente \(MSE\)

## [1] 0.02113112

Finalmente, concluimos que en este caso el boosting no era la mejor opción, pues encontramos que el \(MSE\) para la regresión lineal multiple era mucho menor, y de hecho es el menor de los tres \(MSE\) de prueba calculados.

f) Importancia de las variables

Ahora observaremos cual es la improtancia de las variables respecto a este método. Para ello entrenemos el modelo con el mejor lambda obtenido y analicemos la tabla de importancias.

##                 var    rel.inf
## CAtBat       CAtBat 20.6291342
## Walks         Walks  9.2793701
## PutOuts     PutOuts  8.4533942
## CRBI           CRBI  7.2094278
## CHits         CHits  7.0194877
## CWalks       CWalks  6.8918396
## AtBat         AtBat  5.7101984
## Years         Years  5.1865236
## Assists     Assists  5.1594774
## CHmRun       CHmRun  4.6688414
## HmRun         HmRun  3.9843632
## Hits           Hits  3.7792925
## Runs           Runs  3.6741427
## RBI             RBI  3.0052424
## Errors       Errors  2.1478273
## CRuns         CRuns  2.0991295
## Division   Division  0.6511851
## NewLeague NewLeague  0.3220965
## League       League  0.1290266

Se puede concluir entonces que las variables CAtBat, CHits y CRuns son las variables que más influyentes dentro del modelo.

g) Aplicación de bagging

En este caso aplicaremos bagging sobre este conjunto de datos y observaremos su \(MSE\) de prueba. Dicho \(MSE\) se muestra a continuación

## [1] 0.04434006

Como conclusión podemos observar que el bagging no ofrece un cambio positivo respecto al error, pues al analizar los \(MSE\) de prueba obtenidos con los anteriores métodos, el método lineal sigue siendo el más efectivo de todos.

Sección 9.7 - Ejercicio 4

El objetivo de este ejercicio es verificar la efectividad del metodo SVM para clasificar un conjunto de datos donde no hay una separación lineal entre ellos. Observaremos entonces una comparacion entre un clasificador de soporte vectorial y una maquina de soporte vectorial con kernel polinomial o con un kernel radial.

Primero crearemos el conjunto de datos sintéticos con una separacion no lineal.

## # A tibble: 6 x 3
##      x1    x2 Class
##   <dbl> <dbl> <fct>
## 1 0.587 0.965 0    
## 2 0.139 0.381 0    
## 3 0.301 0.134 0    
## 4 0.866 0.172 0    
## 5 0.622 0.527 1    
## 6 0.394 0.557 1

Grafiquemos el conjunto de datos generado.

Ahora dividiremos el conjunto de datos en un conjunto para entrenamiento y otro conjunto para validación, de la siguiente manera.

## [1] 75  3
## [1] 25  3

Primero entrenemos un SVC (Support Vector Classifier) y extraigamos su error en el conjunto de entrenamiento. Es importante notar que tomaremos el costo para la frontera de clasificación como 10 y observaremos el de este costo sobre los demás modelos.

## 
## Call:
## svm(formula = Class ~ ., data = syntetic_train, kernel = "linear", 
##     cost = 10, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  linear 
##        cost:  10 
## 
## Number of Support Vectors:  72
## 
##  ( 37 35 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

Grafiquemos la clasificación y observemos los resultados.

Y ahora, vamos a extraer el error sobre el conjunto de entrenamiento.

Note entonces que la region clasifica de manera correcta el 53.3333333% de los datos de entrenamiento (esta es la tasa de clasificación correcta de entrenamiento para dicho modelo).

A continuación, computemos la tasa de clasificación correcta para el conjunto de validación.

Se puede notar entonces que para el conjunto de datos de validación, se tiene una tasa de clasificación correcta del 68%. Note que en este caso la grafica es aproximada a la de entrenamiento, pues este clasificador es muy similar a tomar cualquier pareja \(X_1\) y \(X_2\) y asignarles de manera aleatoria uniforme la etiqueta 1 o 0.

Veamos el comportamiento de un kernel radial para este conjunto de datos. Para este caso, tomaremos un valor de \(\gamma = 1\).

## 
## Call:
## svm(formula = Class ~ ., data = syntetic_train, kernel = "radial", 
##     gamma = 1, cost = 10, scale = FALSE)
## 
## 
## Parameters:
##    SVM-Type:  C-classification 
##  SVM-Kernel:  radial 
##        cost:  10 
## 
## Number of Support Vectors:  30
## 
##  ( 16 14 )
## 
## 
## Number of Classes:  2 
## 
## Levels: 
##  0 1

Realicemos una gráfica para este clasificador sobre el conjunto de entrenamiento, como se muestra a continuación.

Es importante apreciar que el clasificador radial hace una clasificación casi perfecta, o al menos así se aprecia graficamente, además toma un conjunto mucho mas razonable de vectores de soporte, que son los que están más cercanos a la frontera de clasificación, al contrario del metodo lineal, el cual tomaba vectores de soporte que no corresponden a la intuición del modelo.

Ahora, para verificar el argumento que se mostró en el parrafo anterior, verifiquemos el error de entrenamiento de este modelo.

Note entonces que la tasa de clasificación correcta para el conjunto de entrenamiento tomando un kernel radial es de 100%, lo cual muestra una efectividad superior al clasificador de soporte vectorial (usando kernel lineal)

De igual manera, observemos la tasa de clasificación correcta para este modelo con kernel radial.

Note entonces que en este caso, la tasa de clasificación correcta es de 92%, lo que muestra que el kernel radial es superior al kernel lineal cuando no hay una separación lineal entre las clases, como era de esperarse. Sin embargo el cambio entre ambos modelos es muy significativo, pues se pasa de tener un modelo cuya selección en la practica corresponde a un clasificador aleatorio uniforme, a tener un modelo con un desempeño casi perfecto.